function measure = measure(f, d)
%MEASURE of a CHEBFUN F betweeen A and B.
%   MEASURE(F, A, B) computes the number F^-1([a,b]) i.e. the measure of
%   the set which is mapped to values between A and B under the mapping F.

% Copyright 2011 by The University of Oxford and The Chebfun Developers. 
% See http://www.maths.ox.ac.uk/chebfun/ for Chebfun information.

% d must be a 2X1 or 1X2 vector
if ( ~isvector(d) || length(d) > 2)
    error('CHEBFUN:MEASURE', 'second argument must be a vector of lenght 2.');
end
a = d(1); 
b = d(2);

% Function f must be real valued:
if ( ~isreal(f) )
    error('CHEBFUN:MEASURE', 'Complex valued functions can not be measured in Chebfun');
end

% If f is not continuous, results may be wrong.
if ( ~iscont(f) )
    warning('CHEBFUN:MEASURE', 'The function is not continuous, results may be inaccurate');
end

% a should not be greater than b:
if ( a > b || a == inf || b == -inf)
    error('CHEBFUN:MEASURE', 'a should not be greater than b' );
end

% What to do if a == b
% if ( a == b ) end

if ( a == -inf && b == inf )
    measure = f.ends(end) - f.ends(1);
    return
end

% If a = -inf and b is bounded:
if ( a == -inf || a < min(f) )
    % A continuous function can not hit a:
    ra = [];
else
    % Compute the points where the function hits a:
    ra = roots(f-a);
end

% If b == inf and a is bounded:
if ( b == inf || b > max(f) )
    % A continuous function can not hit b:
    rb = [];
else
    % Compute the points where the function hits b:
    rb = roots(f-b);
end

% Sort and 'typify' these points into a single array.
[r, t] = sortAndMerge(ra, rb);

% Compute the derivative of f:
fp = diff(f);

% Compute the measure:
measure = 0;
for i = 1:length(r)-1
    if ( t(i) == 'a' )
        if( t(i+1) == 'b' )
            % Function hits a then b
            measure = measure + r(i+1) - r(i);
        else
            % Two consecutive roots hitting a, make sure the function is
            % inside a and b:
            if ( feval(fp, r(i)) >= 0 && feval(fp, r(i+1)) <= 0 )
                measure = measure + r(i+1) - r(i);
            end
        end
    else
        if( t(i+1) == 'a' )
            % Function hits b then a
            measure = measure + r(i+1) - r(i);
        else
            % Two consecutive roots hitting b, make sure the function is
            % inside a and b:
            if ( feval(fp, r(i)) <= 0 && feval(fp, r(i+1)) >= 0 )
                measure = measure + r(i+1) - r(i);
            end
        end
    end
end


% Search for measure at the boundary points.
ends = [f.ends(1), f.ends(end)];
fLeft = feval(f, ends(1));
fRight = feval(f, ends(end));

if ( ~isempty(r) )
    % Check if there is any measure from the left boundary to the first root.
    if ( fLeft >= a && fLeft <= b )
        if ( t(1) == 'a' )
            % Function hits a, check derivative:
            if ( feval(fp, r(1)) <= 0 )
                measure = measure + r(1) - ends(1);
            end
        else
            % Function hits b, check derivative:
            if ( feval(fp, r(1)) >= 0 )
                measure = measure + r(1) - ends(1);
            end
        end
    end
    
    % Check if there is any measure from the last root to the right boundary.
    if ( fRight >= a && fRight <= b )
        if ( t(end) == 'a' )
            % Function hits a, check derivative:
            if ( feval(fp, r(end)) >= 0 )
                measure = measure + ends(end) - r(end);
            end
        else
            % Function hits b, check derivative:
            if ( feval(fp, r(end)) <= 0 )
                measure = measure + ends(end) - r(end);
            end
        end
    end
else
    % Function does not hit a or b.
    if ( fLeft >= a && fLeft <= b )
        if ( fRight >= a && fRight <= b )
            % Measure is the length of the whole interval.
            measure = ends(end) - ends(1);
        else
            % This means the function did hit a or b and but not detected.
            error( 'CHEBFUN:MEASURE', 'root not detected');
        end
    end
end
end

function [r, t] = sortAndMerge(ra, rb)
%SORTANDMERGE constructs a single sorted array obtained from ra and rb.
%   SORTANDMERGE(RA, RB) merges RA and RB and sorts the merged list. Vector
%   T is of the same length as the merged array such that T(i) is either
%   'a' or 'b' indicating the original array from with the ith member of
%   the merged array is coming from. It is assumed that RA and RB have no
%   elements in common.

m = 1; n = 1; i = 1;
lenA = length(ra);
lenB = length(rb);

r = zeros(lenA + lenB, 1);
t = zeros(lenA + lenB, 1);

while ( m <= lenA || n <= lenB )
    % If one of the arrays is exausted, copy the second one into output:
    if ( m > lenA )
        while ( n <= lenB )
            r(i) = rb(n);
            t(i) = 'b';
            n = n + 1;
            i = i + 1;
        end
        break;
    end
    
    if ( n > lenB )
        while ( m <= lenA )
            r(i) = ra(m);
            t(i) = 'a';
            m = m + 1;
            i = i + 1;
        end
        break;
    end
    
    % None of the arrays is exausted;
    if ( ra(m) < rb(n) )
        r(i) = ra(m);
        t(i) = 'a';
        i = i + 1;
        m = m + 1;
    else
        r(i) = rb(n);
        t(i) = 'b';
        i = i + 1;
        n = n + 1;
    end
end
% Sanity check
if ( length(r) ~= lenA + lenB )
    % Something has gone wrong.
    error('CHEBFUN:HIST:SORTANDMERGE', 'Roots likely to have duplication');
end

end